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ABSTRACT 

Recent studies have indicated that large mass radial inflow rates are possible in ob- 
served galaxies with strong density wave patterns. Yet, numerical simulations have 
generally failed to account for such high rates. Here it is shown that the reason for 
the discrepancy is the treatment of "softening", the artificial parameter inserted by 
numerical simulators into the formula for gravitational potential, to control the mag- 
nitude of relaxation in simulations with small numbers of particles compared to a real 
galaxy. Excess softening reduces the collective effects underlying the significant secu- 
lar evolution inferred for physical galaxies. Less softening, coupled with an increase in 
number of particles, allow the N-body simulations to reveal significant morphological 
transformation of galaxies over a Hubble time. 
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1 INTRODUCTION 

The prospect of significant transformation of the Hubble 
type of a galaxy during its lifetime has important impli- 
cations on our understanding of the origins of the Hubble 
sequence and on broader issues in cosmology. Despite the 
coming of age of this field of research, there is as of now 
no consensus on how important this process is, and through 
what dynamical mechanism this transformation is mainly 
achieved (see the review of Kormendy & Kennicutt 2004) . 

Zhang (1996, 1998, 1999) proposed and demonstrated 
a collective dissipation process in galaxies that could serve 
as the underlying dynamical mechanism for secular evolu- 
tion of galaxies along the entire Hubble sequence, by incor- 
porating the stellar component in the mass inflow process 
rather than just the interstellar medium as proposed in the 
pseudo-bulge formation scenario (Kormendy 1979). Despite 
the confirmation of the theoretically derived mass flow rates 
in the accompanying N-body simulations, which established 
the viability of the analytical approach, the simulated mass 
flow rates were small, far from being adequate to lead to 
significant morphological transformation of a galaxy within 
a Hubble time. This small mass flow rate was attributed to 
the small amplitudes of the spiral patterns formed in these 
simulations, with the state-of-the-art in simulating sponta- 
neously formed modes far from being adequate to produce 
the extremely nonlinear density wave modes observed in 
physical galaxies. 

In order to circumvent this difficulty, Zhang & Buta 



(2007; 2012) used near-and-mid-infrared images of galaxies, 
coupled with the analytical rate equations derived in Zhang 
(1996, 1998, 1999) to calculate the torques and mass flow 
rates in these physical galaxies, and found that mass flow 
rates from a few Mq per year to over one hundred Mq per 
year were obtained, depending on the Hubble types of galax- 
ies and their interacting environment, which together deter- 
mined the density wave amplitudes and pitch angles in these 
galaxies. These mass flow rates for real galaxies far exceed 
those obtained from the past N-body simulations, partly as 
a result of the fact that the mass flow rate is proportional to 
the effective wave amplitude squared (Zhang 1998), and the 
amplitudes in physical galaxies are oftentimes a factor of 3-8 
times higher than obtained in N-body simulations, implying 
a difference in evolution rate of a factor of 10 - 100 between 
physical and previously-simulated galaxies. 

Despite the confirmation of the importance of secular 
evolution process in physical galaxies in works such as Zhang 
& Buta (2007,2012), one could not help but wonder what 
exactly were the factors that have prohibited the simulated 
disk galaxies from obtaining the level of mass flow rates and 
the magnitude of wave amplitudes in physical galaxies. It is 
to this question that the result of the current Letter partly 
addresses. Incidentally, the main discovery of this work was 
made in a totally serendipitous fashion - i.e. the result was 
stumbled upon while the author was searching in a different 
direction for the cures of the low mass flow rates in N-body 
simulation^]) . The lesson learned is that one has to be cau- 
tious in taking the results of N-body simulations literally 
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1 The author used to think that the low mass flow rates were 
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without questioning sacred "rules of thumb" and accepted 
wisdom, especially in circumstances where new physical pro- 
cesses are encountered - in our situation, that is when the 
collective effects become important. 



2 ROLES OF SOFTENING: THEORETICAL 
AND NUMERICAL 

To come to the core result of this Letter immediately: the 
main culprit to the low mass accretion rate in N-body simu- 
lations of disk galaxies turned out to be the improper choice 
of softening parameter used to control relaxation effects in 
the simulation of these nominally collisionless systems. 

The softening parameter is generally defined as follows: 
The gravitational potential in an N-body system is calcu- 
lated as 



where a so ft is the so-called softening parameter. This pa- 
rameter is introduced into the N-body simulations due to 
the fact that the numerical setups usually have significant 
less number of particles than the physical systems they try to 
emulate. A finite softening parameter reduces the amount of 
unrealistic close encounters which would be more often for a 
small N system, and thus reduces the level of artificial relax- 
ation in these small N systems. The hope is that a collision- 
less configuration would then result from a well-matched N 
and a so ft. In particle- mesh based N-body approaches, the 
grid size serves as an additional softening parameter that 
acts in consort with a so ft in controlling relaxation (see Sell- 
wood 1987 and the references therein). 

However, with the realization that galaxies which pos- 
sess large scale density wave modes are governed by col- 
lective dissipation processes (Zhang 1996), one soon faces 
the fact that it is impossible to completely avoid collision- 
like behavior, no matter how large the particle numbers are 
used: After all, the collective effects in these systems depend 
on the near-collision or scattering of particles in the global 
instabilities to set up long-range correlations to achieve self- 
organization and to induce secular evolution. In some sense, 
these galactic systems are in a forced relaxation configura- 
tion, with the forcing accomplished by the density wave col- 
lisionless shocks, which force the effective Q in the spiral 
arms to be less than one and lead to inter-particle interac- 
tions and correlations. Thus, in hindsight, it is also obvious 



mainly due to the two-dimensional (2D) nature of her simulations, 
and once the bulge and halo were allowed to be active, the mass 
flow rates will become higher. This turned out to be true, but only 
moderately true: i.e., having about 15% active bulge and halo in 
2D simulations increases the mass inflow rate by about 10%, but 
having still larger active bulge and halo fraction leads rather to 
outflows of mass over the main parts of the disk (likely due to the 
lowering of Q and the emergence of transient local instabilities), 
contrary to what is needed for bulge building. Note that the above 
conclusion about the fraction of allowable active bulge/halo mass 
applies only in the 2D cases where these tests had been made. 
In three dimensional (3D) studies the entire bulge and halo mass 
of course can be made active with proper pressure and rotational 
support. 



that excess softening reduces the very inter-particle inter- 
action that is the backbone support of collective effects. In 
essence, overly-softened gravity is modified gravity, and thus 
is modified Newtonian interaction (not to be confused with 
the Modified Newtonian Dynamics proposed to account for 
the flat rotation curves) . The system one simulates with soft- 
ened gravity is no longer exactly the same system governed 
by Newton's gravitational law, but rather a more sluggishly- 
interacting system. This might not be a serious concern if 
one's interest is in obtaining a modal morphology that mim- 
ics the observed galaxy morphology (as shown in the work 
of Donner & Thomasson 1994, hereafter DT94, which is the 
forerunner of long- lasting N-body spiral modes), but as we 
will show here, it is detrimental to determining the realistic 
secular mass flow rates that are relevant to physical galaxies. 

In Fig. [TJ we show the result of a set of test runs of 
N-body simulations using the basic state first explored in 
DT94, and subsequently used for the study of collective ef- 
fects and secular evolution in Zhang (1996,1998,1999). The 
simulations are done in a 2D configuration, using a particle- 
mesh approach on a polar grid, with an active disk which 
has a modified exponential surface density profile, as well as 
a rigid bulge and a rigid halo. The grid size is 220 by 256 in 
the radial and azimuthal direction, respectively. The mass 
in the active disk, the inert halo and inert bulge are 0.4, 0.5, 
0.1, respectively, in the normalized unit, similar to that used 
in Zhang (1998). The scaling of the time unit is such that 
1256 time steps represent one rotation period at r=20. Both 
the particle numbers and the time resolution used are higher 
than that in Zhang (1998). These grid and time step reso- 
lutions are used for the rest of the simulations presented in 
this Letter as well, only particle numbers and the softening 
paramters are changed and will be noted in each case. 

The curves in Fig.l are for the evolution of enclosed 
mass within the central r — 3.5 region (in the scale of this 
disk which has corotation radius around 20-25 in normalized 
unit, this corresponds to the central bulge region). The three 
curves, from bottom to top, correspond to softening parame- 
ter value of 1.5 (as used in DT94 and Zhang 1998), 0.75, and 
0.25 in the unit of the smallest grid size of 1. As we can see, 
the choice of softening of 1.5 produced a barely noticeable 
mass inflow (to really see the small amount of mass inflow for 
this choice of softening, one is referred to Figure 10 of Zhang 
1998 which zoomed in on the vertical scale, even though that 
was for the enclosed mass within central r = 13.5). For that 
particular choice of softening, which is of standard usage 
in the N-body simulation community as being adequate in 
suppressing unwanted relaxation effects as well as in match- 
ing the grid softening, the central mass growth of only 6.5% 
was observed over about 25 pattern rotation periods (Zhang 
1998), which is obviously quite inadequate in transforming 
the Hubble type of a galaxy. Gradually reducing softening is 
seen to systematically lift up the mass inflow rate. For the 
range of softening tried by the current author (1.5-0.1), the 
trend of increasing mass inflow with reduced softening con- 
tinues. It is possible that realistic galaxies, which have zero 
softening in the form of gravitational law but do have much 
larger numbers of particles, as well as the finite thickness 
of the disk to keep relaxation in check, the mass flow rates 
much exceed those we have seen in these recent simulations 
are possible - as we indeed observed in the mass flow rates 
calculated for the sample in Zhang & Buta (2007, 2012). 



© 0000 RAS, MNRAS 000, 000-000 



Rapid Secular Evolution in N-Body Galaxy Disks 3 




o ~_ 
X 

CO - 

I 

O " 
X 

S' X - 
CO - 

I 

o ~_ 

X - 
m - 

CT3 



5000 10 4 1.5xl0 4 2xl0 4 

Time Step 

Figure 1. Evolution of enclosed disk mass within central r=3.5 
radius of a set of 2D N-body simulations with different softening 
parameters. Solid: a so f t = 1.5, Dashed: a so f t = 0.75, Dotted:i 
&soft = 0.25. 1 million particles are used to represent the active 
disk. 



3 A NEW SET OF N-BODY SIMULATIONS 
WITH SMALL SOFTENING: 
ACCELERATED BULGE BUILDING 

The price one pays for reduced softening in small-N simula- 
tions is the increase in relaxation rate. To compensate, one 
needs to use a correspondingly increased number of parti- 
cles, for otherwise the collision-induced heating will swamp 
the density wave activity and reduce the wave amplitude. 
How much the particle number needs to be increased for a 
given choice of softening can be found empirically (recalling 
the fact that the traditional rules of thumb in cases with col- 
lective effects may not always apply, since we do not want to 
kill the right kind of collisional relaxation due to collective 
effects). 

In what follows, we show a set of simulations with 
changing particle numbers while holding the softening pa- 
rameter a so f t =0.1. Fig. [2] shows the morphological evolu- 
tion of a spontaneously formed spiral-bar mode for one of 
these runs using 20 million particles. Note that even though 
the duration of run for this mode is similar to that used in 
Zhang (1998) (note due to the different temporal resolution 
the 400 step here corresponds to 100 step in the older simu- 
lation) , the modal morphology here appears to have changed 
much more that that seen in Zhang (1998), i.e. the pattern in 
the older simulation retained the spiral morphology through- 
out, whereas the pattern in the current simulation evolved 
from a spiral- like morphology to a more bar- like morphol- 
ogy. This is natural to expect, since the much-higher radial 
mass accretion rate in the current simulation resulted in a 
significant increase in central mass concentration at the end 
of the run, this further means that the basic state is change 
so much that the mode it supports needs also to change 
correspondingly from that of spiral to that of bar with pos- 
sibly nested inner spiral (note that a heavy bar is the nat- 
ural modal shape for a centrally- heavy disk - the observed 
early- types disks certainly contain many massive bars) . The 
pattern speed, however, appears to be little affected by the 
change in morphology, as can be seen later in the clean con- 
centration of the power spectrum around a single pattern 
speed throughout the duration of the run. 

In Fig. [3l we present the enclosed mass within the cen- 
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Figure 2. N-body morphology of a spiral/bar mode at different 
time steps. The rotation period at r=20 is about 1256 steps. 20 
million particles are used to represent the active disk, the soft- 
ening parameter is a so f t =0.1. About 10,000 randomly-selected 
particles are plotted. 



tral r=3.5 radius, with the different runs having a constant 
a so ft = 0.1 but varying the number of particle used. For the 
smaller particle number runs, the rapid mass inflow is seen 
to saturate at an earlier time step, which corresponds to the 
time when spiral activity is damped by heating due to the 
insufficient number of particles for this extremely small soft- 
ening, especially in the outer disk region where the surface 
density is low and the particle numbers are small. For large 
particle runs, the rapid mass inflow is seen to remain at a 
constant rate all the way until the end of the run, which 
corresponds to 25 rotations. 

Furthermore, we note that the mass inflow rate observed 
here corresponds to about 60% increase in enclosed mass 
within r=3.5 (the bulge region) over 10 rotation periods, 
or roughly 1/5 of a Hubble type. This level of mass accre- 
tion is more than sufficient to transform the Hubble type 
of a galaxy by several stages in a Hubble time. Plus, the 
mass inflow rate conceivably could be even higher for even 
smaller choices of a so ft- The natural limit apparently is set 
only when one has the actually physical galaxy in its 3D 
representation with the full inter-particle interactions taken 
into account (i.e. we know what the actual density wave 
amplitudes in galaxies are because we observed them). 

In Fig. [4] the disk surface density at the beginning of 
the run, and at the end of the four runs with different num- 
bers of particles (as used in Fig [3] but with different line 
style due to the need to represent the "before" surface den- 
sity) are presented. It can be seen that the large inflow of 
the large-particle-number runs indeed builds up a more sub- 
stantial bulge. However, the 20 million particle run in fact 
produced more substantial bulge-building than the 40 mil- 
lion particle run. This is consistent with the trend observed 
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Figure 3. Enclosed mass within central r=3.5 of the N-body 
disk for runs with different numbers of particles and a s oft = 0.1. 
Solid: 1 million particles. Dashed: 10 million particles. Dotted: 20 
million particles. Dash-Dotted: 40 million particles. 
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Figure 4. Basic state surface densities at the beginning (solid), 
and the end of four N-Body runs with different number of Par- 
ticles. Dashed: end surface density for 1 million particles. Dot- 
ted: end surface density for 10 million particles. Dash-and-Single- 
Dotted: end surface density for 20 million particles. Dash-and- 
Triple-Dotted: end surface density for 40 million particles. a so f t = 
0.1 

in Fig. [3l where it is seen that the 40 million particle run 
has somewhat reduced mass inflow rate near the end of the 
run. This may be related to the emergence of different kinds 
of local instabilities in the different particle runs (a small 
knot of globular- cluster- like morphology can be seen in the 
outer disk of the second and third frame of Fig. [2]). This 
kind of cluster-forming activity in association with density 
wave modal activity may be a reflection of the actual process 
happening in physical galaxies. By step 32000 (not shown in 
Fig [2]), the 20 million particle run has developed two promi- 
nent dumb-bell like structures at the end of the massive bar 
which might account for the increased mass inflow for this 
run. 

In Fig. [5] the m=2 power spectrum together with the 
Q, Q±k/2 curves are presented for two runs with two differ- 
ent numbers of particles. These contour plots clearly show 
a single dominant mode throughout the time interval of the 
runs. These are the cleanest power spectra (albeit obtained 
using potential m=2 components rather than the more con- 




Figure 5. Power spectra from m=2 Potential at different radii 
during time interval step=0 and step=32768. Top: 20 million par- 
ticle run. Bottom: 40 million particle run. The solid line indicates 
galactic rotation speed fl, the two dashed lines are fl ± k/2. 



ventional surface density) yet observed for this particular 
mode, compared to those runs with larger softening and 
smaller number of particles (note that due to higher tem- 
poral resolution by a factor of 4 compared to that used in 
Zhang 1998, the power spectrum should be multiplied by 
by the same factor of 4 to compare with older results. Af- 
ter taking this into account, the pattern speed still appears 
to be slightly higher than before - even the 40 million run 
appears to have a slightly higher pattern speed than the 20 
million particle run). 

In Fig. [6] we plot the mass inflow (positive part of the 
curve) and outflow (negative part of the curve) versus ra- 
dius at time step 8000. This mass flow rate curve was cal- 
culated based on the volume torque prescription given in 
Zhang (1996, 1998). The major transition from inflow to 
outflow is close to r=20 for the inner disk, which is close to 
(but not exactly the same, due probably to the fast evolving 
nature of the surface density) the corotation radius predicted 
from pattern speed (Fig. [5}. Note that the larger number-of- 
particle case the power spectrum predicts a slightly smaller 
corotation radius than the smaller number-of-particle run. 
The level of peak accretion, 2 — 4 x 10 -7 per time step for 
the inner disk region, matches quite well the average ac- 
cretion rate shown from the slope of the Fig. [3] 20 million 
particle curve, which corresponds to 2 x 10 -7 per time step, 
within the central r=3.5. It should be noted that the small 
peaks and valleys of the mass flow curve in the inner disk 
do change with time and in general the peaks move inward 
with time. 

Thus this new set of simulations show that the broad 
modal nature of the pattern as well as the predictions of 
the collective dissipation formalism persisted even with the 
drastic reduction in softening length. The dominant mode 
morphology does evolve a lot faster in appearance, which is 
to be expected in a fast changing basic state. 

Finally we note that with the drastically-reduced a so ft 
in the equation for the N-body gravitational potential in our 
new set of simulations, the corresponding grid softening has 
not been much reduced (the newer grid is only about a fac- 
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Figure 6. Mass inflow (positive) and outflow (negative) versus 
radius at time step 8000 

tor of 4 finer in linear resolution than used in Zhang 1998, 
whereas the softening is reduced by a factor of 15). In fact 
Zhang 1996 already explored a grid twice as fine as used in 
DT94 (or in Zhang 1998) and the mass inflow rate was not 
found to have changed much. Thus it appears that particle 
softening is the main inhibitor to obtaining higher mass in- 
flow rates in N-body disks containing collective modes, the 
grid softening has only marginal effect in this regard. This 
is reasonable as the main contributor to collective effects 
is inter-particle interactions, and the finite amount of grid 
softening perhaps also helped to maintain the desirable finite 
thickness effect in these 2D simulations. 



4 DISCUSSION 

4.1 Are the Increased Mass Flow Rates in 
Less-Softened N-Body Runs Due to the 
Collective Instabilities or Due to Few-Body 
Encounters? 

The answer is that it is due to collective instabilities asso- 
ciated with the density wave modes. The hint for this con- 
clusion can be found again in Fig. [3] It is seen there that 
when not enough number of particles are used for a given 
(isoft (which generally increases the few-particle relaxation 
effect), the mass inflow tapers off at an earlier time step. 
Whereas for very large particle numbers (which generally 
decreases the few-particle relaxation effect), both the spi- 
ral/bar activity and the mass inflow continues until the end 
of the run. This correlation between the number of simu- 
lation particles and the longevity of effective mass inflow 
shows that the survivability of the density wave modes is 
key to the continued mass inflow. Few-particle effect causes 
negligible mass inflow once the spiral/bar activity ceases. 

4.2 The Remaining Unrealistic Elements in the 
Secular Bulge Building Simulations 

The simulations described in this Letter serve to highlight 
an important, and so far neglected, aspect of N-body simu- 
lations that has a huge impact on the ability of these sim- 
ulations to reproduce the observed level of mass inflow for 
bulge building as inferred from the MIR observations (Zhang 



& Buta 2007, 2012). However, since the simulation is con- 
ducted in 2D, obvious physical effects that are missed in- 
clude the continued growth of the disk through halo mass 
accretion, the role of active disk and halo (though as dis- 
cussed in the footnote, we have a general idea of what this 
will do to the enhancement of the mass inflow rate). Fur- 
thermore, since this is a pure particle simulation, the role 
of gas is ignored, though the star-gas two component sim- 
ulation conducted in Zhang (1998) shows that the general 
features of the collective effects are preserved when gas is 
added, and Zhang & Buta (2012) showed that the overall 
effect of stellar accretion is more important than gas accre- 
tion due to the larger stellar surface density, for all galaxies 
except the very late types. 



5 CONCLUSIONS 

The low mass inflow rates found in the past N-body simula- 
tions of disk galaxies had cast doubt on the effectiveness of 
secular evolution as an important dynamical process for the 
morphological transformation of galaxies along the Hubble 
sequence. In this work it is shown that such low numerical 
mass inflow rates are chiefly the result of the artificial "soft- 
ening" of gravity which is a common practice in galaxy sim- 
ulations to avoid rapid relaxation due to the small number 
of particles employed. By decreasing the amount of soften- 
ing and simultaneously increasing the number of simulation 
particles, realistic level of mass inflow comparable to those 
inferred from observed galaxy images through torque calcu- 
lations are achieved in N-body simulations of disk galaxies, 
which reinforce the importance of secular evolution as an ex- 
tremely relevant process in transforming the morphologies 
of galaxies during the past Hubble time. 
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